using System;
using System.Collections.Generic;
using System.Text;
using interfaces;

namespace retrofit
{
    public class Algorithm
    {
        private static int m_bilinearNum = 2;

//		private static double m_variationConst = 0.2;	// variation constant is defined for tables to compare with baseline case.

		private Query m_query;

		enum Refrigeration_T
		{
			Refrigeration_RT,		// no refrigeration case, HasRefrigeration = 0
			Baseline_RT,			// baseline case with refrigeration, HasRefrigeration = 1.
			NUM_RT
		};


        /// <summary>
        /// constructor
        /// </summary>
        public Algorithm(ref Query q)
        {
			m_query = q;
        }

		
        //** ---------- process section ------------ **//
        
		/// <summary>
		/// process all energy types.
		/// </summary>
		/// <param name="tables"> baseline energy output array for all table types<param name="Output"></param></param>
		/// <return> reference of calculation output</return>
		public Output process(string building_type, Input input, double[] e_factors, double[] global, List<string> err_list)
		{
			List<Output> outputs = new List<Output>();
			for (int i = 0; i < (int)Keyword.table_t.NUM_TT; ++i)
			{
				Output out1 = new Output();
				outputs.Add(out1);
			}

			processDiscreteTables(building_type, input, outputs, err_list);

			processContinuousTables(building_type, input, outputs, err_list);

//			processEnvelopeTable(input, outputs, err_list);

			Output out2 = new Output();

			List<double[]> cross_factors = refactorTables(outputs);

			// process HVAC energy
			processHVAC(outputs, e_factors, cross_factors, out2);

			// process other energy
			processOthers(outputs, e_factors, out2);

//			Console.WriteLine("BEFORE");
//			printArray(out2.m_energy);

			Output out_final = processGlobalFactors(out2, global);

//			Console.WriteLine("AFTER:");
//			printArray(out_final.m_energy);

			return out_final;
		}


        /// <summary>
        /// process other hvac related energy.
        /// </summary>
        /// <param name="tables"> baseline energy output array for all table types<param name="Output"></param></param>
		/// <param name="out1"> reference of calculation output</param>
        private static void processHVAC(List<Output> tables, double[] e_factor, List<double[]> cross_factor, Output out1)
        {
            //** ---- first we process all heat load table (envelope, roof, int lighting, int equipment) for HVAC energy ---- **//

			// TODO: error handling

            // energy difference, and each difference's scale
            Output diff, scale;  
            diff = new Output();
            scale = new Output();
            
            // we only loop those tables that affect building's internal heat load.
			// TODO, consider lighting and equipment internal load. Integrate with empirical
			for (int i = 0; i <= (int)Keyword.table_t.Refrigeration_TT; ++i)
            {
				// we only loop through all the HVAC components
				for (int j = 0; j <= (int)Keyword.energy_cat_t.CoolingTower_ECT; ++j)
				{
					// for internal heat load energy, which can affect HVAC energy, we use cross_factro to scale it.
					if (i >= (int)Keyword.table_t.IntLight_TT)
					{
						diff.m_energy[j] -= tables[(int)Keyword.table_t.Baseline_TT].m_energy[j] * e_factor[(int)Keyword.m_factorBaseEnergy[i]] * cross_factor[i][j];
					}
					// for envelope and other table which is query based results, we just use what we get from query.
					else
					{
						diff.m_energy[j] += tables[i].m_energy[j] - tables[(int)Keyword.table_t.Baseline_TT].m_energy[j];
					}
				}
			}

			// we scale each HVAC energy based on baseline number
			for (int i = 0; i <= (int)Keyword.energy_cat_t.CoolingTower_ECT; ++i)
			{
				scale.m_energy[i] = diff.m_energy[i] / tables[(int)Keyword.table_t.Baseline_TT].m_energy[i];
				//Console.WriteLine(diff.m_energy[i] + " ," + scale.m_energy[i]);
			}

			//** ---- then, we use the scale to process HVAC energy ---- **//
			for (int i = 0; i <= (int)Keyword.energy_cat_t.CoolingTower_ECT; ++i)
			{
				// process heating
				if (i == (int)Keyword.energy_cat_t.Heating_ECT)
					out1.m_energy[i] = tables[(int)Keyword.table_t.Heating_TT].m_energy[i] * (1 + scale.m_energy[i]) * (1 - e_factor[i]);
				// process fan
				else if (i == (int)Keyword.energy_cat_t.Fan_ECT)
					out1.m_energy[i] = tables[(int)Keyword.table_t.Fan_TT].m_energy[i] * (1 + scale.m_energy[i]) * (1 - e_factor[i]);
				// process cooling, pump and CT
				else
					out1.m_energy[i] = tables[(int)Keyword.table_t.Cooling_TT].m_energy[i] * (1 + scale.m_energy[i]) * (1 - e_factor[i]);
            }
        }


		/// <summary>
		/// process other energy formate other than HVAC.
		/// </summary>
		/// <param name="tables"> baseline energy output array for all table types<param name="Output"></param></param>
		/// <param name="out1"> reference of calculation output</param>
		private static void processOthers(List<Output> tables, double[] e_factor, Output out1)
        {
            // refrigeration
			out1.m_energy[(int)Keyword.energy_cat_t.Refrigeration_ECT] = (1 - e_factor[(int)Keyword.energy_cat_t.Refrigeration_ECT]) * 
				tables[(int)Keyword.table_t.Baseline_TT].m_energy[(int)Keyword.energy_cat_t.Refrigeration_ECT];

            // interior lighting
			out1.m_energy[(int)Keyword.energy_cat_t.IntLighting_ECT] = (1 - e_factor[(int)Keyword.energy_cat_t.IntLighting_ECT]) *
				tables[(int)Keyword.table_t.Baseline_TT].m_energy[(int)Keyword.energy_cat_t.IntLighting_ECT];

            // exterior lighting
			out1.m_energy[(int)Keyword.energy_cat_t.ExtLighting_ECT] = (1 - e_factor[(int)Keyword.energy_cat_t.ExtLighting_ECT]) *
				tables[(int)Keyword.table_t.Baseline_TT].m_energy[(int)Keyword.energy_cat_t.ExtLighting_ECT];

            // interior equipment ele
			out1.m_energy[(int)Keyword.energy_cat_t.IntEquipmentEle_ECT] = (1 - e_factor[(int)Keyword.energy_cat_t.IntEquipmentEle_ECT]) *
				tables[(int)Keyword.table_t.Baseline_TT].m_energy[(int)Keyword.energy_cat_t.IntEquipmentEle_ECT];

            // interior equipment gas
			out1.m_energy[(int)Keyword.energy_cat_t.IntEquipmentGas_ECT] = (1 - e_factor[(int)Keyword.energy_cat_t.IntEquipmentGas_ECT]) *
				tables[(int)Keyword.table_t.Baseline_TT].m_energy[(int)Keyword.energy_cat_t.IntEquipmentGas_ECT];

			// transportation
			out1.m_energy[(int)Keyword.energy_cat_t.Transportation_ECT] = (1 - e_factor[(int)Keyword.energy_cat_t.Transportation_ECT]) * 
				tables[(int)Keyword.table_t.Baseline_TT].m_energy[(int)Keyword.energy_cat_t.Transportation_ECT];
        }


		public Output processGlobalFactors(Output out_in, double[] global)
		{
			// TODO: validate each global factor is between the range of [-1, 1]

			Output output = new Output();

			// interior lighting
			output.m_energy[(int)Keyword.energy_cat_t.IntLighting_ECT] = (1 + global[(int)interf.BEE_EnergyFactor_T.Lighting_EFT]) *
				out_in.m_energy[(int)Keyword.energy_cat_t.IntLighting_ECT];

			// exterior lighting
			output.m_energy[(int)Keyword.energy_cat_t.ExtLighting_ECT] = (1 + global[(int)interf.BEE_EnergyFactor_T.Lighting_EFT]) *
				out_in.m_energy[(int)Keyword.energy_cat_t.ExtLighting_ECT];

			// interior equipment ele
			output.m_energy[(int)Keyword.energy_cat_t.IntEquipmentEle_ECT] = (1 + global[(int)interf.BEE_EnergyFactor_T.IntEquip_EFT]) *
				out_in.m_energy[(int)Keyword.energy_cat_t.IntEquipmentEle_ECT];

			// interior equipment gas
			output.m_energy[(int)Keyword.energy_cat_t.IntEquipmentGas_ECT] = (1 + global[(int)interf.BEE_EnergyFactor_T.IntEquip_EFT]) *
				out_in.m_energy[(int)Keyword.energy_cat_t.IntEquipmentGas_ECT];

			// transportation
			output.m_energy[(int)Keyword.energy_cat_t.Transportation_ECT] = (1 + global[(int)interf.BEE_EnergyFactor_T.IntEquip_EFT]) *
				out_in.m_energy[(int)Keyword.energy_cat_t.Transportation_ECT];

			// DHM
			output.m_energy[(int)Keyword.energy_cat_t.DHW_ECT] = (1 + global[(int)interf.BEE_EnergyFactor_T.DHW_EFT]) *
				out_in.m_energy[(int)Keyword.energy_cat_t.DHW_ECT];

			// cooling
			output.m_energy[(int)Keyword.energy_cat_t.Cooling_ECT] = (1 + global[(int)interf.BEE_EnergyFactor_T.AirCond_EFT]) *
				out_in.m_energy[(int)Keyword.energy_cat_t.Cooling_ECT];

			// cooling tower
			output.m_energy[(int)Keyword.energy_cat_t.CoolingTower_ECT] = (1 + global[(int)interf.BEE_EnergyFactor_T.AirCond_EFT]) *
				out_in.m_energy[(int)Keyword.energy_cat_t.CoolingTower_ECT];

			// pumps
			output.m_energy[(int)Keyword.energy_cat_t.Pump_ECT] = (1 + global[(int)interf.BEE_EnergyFactor_T.AirCond_EFT]) *
				out_in.m_energy[(int)Keyword.energy_cat_t.Pump_ECT];

			// fans
			output.m_energy[(int)Keyword.energy_cat_t.Fan_ECT] = (1 + global[(int)interf.BEE_EnergyFactor_T.AirCond_EFT]) *
				out_in.m_energy[(int)Keyword.energy_cat_t.Fan_ECT];

			// heating
			output.m_energy[(int)Keyword.energy_cat_t.Heating_ECT] = (1 + global[(int)interf.BEE_EnergyFactor_T.Heating_EFT]) *
				out_in.m_energy[(int)Keyword.energy_cat_t.Heating_ECT];

			// refrigeration
			output.m_energy[(int)Keyword.energy_cat_t.Refrigeration_ECT] = (1 + global[(int)interf.BEE_EnergyFactor_T.Refrig_EFT]) *
				out_in.m_energy[(int)Keyword.energy_cat_t.Refrigeration_ECT];


			return output;
		}



		// TODO: debug this function for all the tables
		private static List<double[]> refactorTables(List<Output> tables)
		{ 
			List<double[]> factors = new List<double[]>();

			for (int i = 0; i < (int)Keyword.table_t.Baseline_TT; ++i)
			{
				double[] f = new double[(int)Keyword.energy_cat_t.NUM_ECT];
				
				// if this table needs to factor, factor each energy category based on its based energy.
				if (Keyword.m_factorBaseEnergy[i] != Keyword.energy_cat_t.NUM_ECT)
				{
					// each table's base energy difference compared with baseline table, we need percentage value of base energy
					double base_energy = tables[(int)Keyword.table_t.Baseline_TT].m_energy[(int)Keyword.m_factorBaseEnergy[i]];
					double table_energy = tables[i].m_energy[(int)Keyword.m_factorBaseEnergy[i]];

					double diff = (tables[i].m_energy[(int)Keyword.m_factorBaseEnergy[i]] - tables[(int)Keyword.table_t.Baseline_TT].m_energy[(int)Keyword.m_factorBaseEnergy[i]])
									/ tables[(int)Keyword.table_t.Baseline_TT].m_energy[(int)Keyword.m_factorBaseEnergy[i]];

					for (int j = 0; j < (int)Keyword.energy_cat_t.NUM_ECT; ++j)
					{
						// how base energy change affects each energy category.
						// the results is: %(energy change) / %(base energy change)
						double change = (tables[i].m_energy[j] - tables[(int)Keyword.table_t.Baseline_TT].m_energy[j]) / tables[(int)Keyword.table_t.Baseline_TT].m_energy[j];
						f[j] = change / diff;
					}
					factors.Add(f);
				}
				// if current table does not need to factor
				else
					factors.Add(f);
			}

			return factors;
		}


		/// <summary>
		/// process refrigeration when refrigeration is not exist. Call this function after calling process() to modify energy based on refrigeration condition.
		/// </summary>
		/// <param name="out1">the processed energy output results. The output is get by calling process() function</param>
		/// <param name="err_list">error message list</param>
		public void processRefrigeration(string building_type, Output out1, List<string> err_list)
		{
			// baseline non-refrigration table
			int[] idx = new int[] { (int)Refrigeration_T.Refrigeration_RT };
			Output non_refrig_out = m_query.queryTable(Keyword.table_t.Baseline_TT, building_type, idx, err_list);

			// baseline table
			idx = new int[] { (int)Refrigeration_T.Baseline_RT };
			Output baseline_out = m_query.queryTable(Keyword.table_t.Baseline_TT, building_type, idx, err_list);


			// we scale each HVAC energy based on baseline (w/ & w/o refrigeration) number
			double scale;
            for (int i = 0; i <= (int)Keyword.energy_cat_t.CoolingTower_ECT; ++i)
            {
				// we only compare heating energy and pump energy on baseline and non-refrigeration baseline.
				if (i == (int)Keyword.energy_cat_t.Heating_ECT || i == (int)Keyword.energy_cat_t.Pump_ECT)
				{
					scale = non_refrig_out.m_energy[i] / baseline_out.m_energy[i];
					out1.m_energy[i] *= scale;
				}
            }

			// clean up refrigeration.
			out1.m_energy[(int)Keyword.energy_cat_t.Refrigeration_ECT] = 0;
		}


        /// <summary>
        /// process discrete table which does not relate to continuous variable input (intlight, extlight, refrigeration ...)
        /// </summary>
        /// <param name="input">input object</param>
        /// <param name="out_list">energy output list coitaning query energy data from variaous tables</param>
        /// <param name="err_list">error list message</param>
		private void processDiscreteTables(string building_type, Input input, List<Output> out_list, List<string> err_list)
		{
			Output out_temp = new Output();
			Output out1 = new Output();

			// exterior lighting table.
			out_list[(int)Keyword.table_t.ExtLighting_TT] = m_query.queryTable(Keyword.table_t.ExtLighting_TT, building_type, input.m_inputIndex, err_list);
			
			// interior equipment table.
			out_list[(int)Keyword.table_t.IntEquipment_TT] = m_query.queryTable(Keyword.table_t.IntEquipment_TT, building_type, input.m_inputIndex, err_list);
			
			// interior lighting table.
			out_list[(int)Keyword.table_t.IntLight_TT] = m_query.queryTable(Keyword.table_t.IntLight_TT, building_type, input.m_inputIndex, err_list);
			
			// refrigeration table.
			out_list[(int)Keyword.table_t.Refrigeration_TT] = m_query.queryTable(Keyword.table_t.Refrigeration_TT, building_type, input.m_inputIndex, err_list);

			// fan table
			out_list[(int)Keyword.table_t.Fan_TT] = m_query.queryTable(Keyword.table_t.Fan_TT, building_type, input.m_inputIndex, err_list);

			// baseline table, to distingush with/without refrigeration case for shoppingmall building
			int[] idx = new int[] { (int)Refrigeration_T.Baseline_RT };
			out_list[(int)Keyword.table_t.Baseline_TT] = m_query.queryTable(Keyword.table_t.Baseline_TT, building_type, idx, err_list);
		}


        /// <summary>
        /// process continuous table except building envelope table
        /// </summary>
        /// <param name="input">input object</param>
        /// <param name="out_list">energy output list coitaning query energy data from variaous tables</param>
        /// <param name="err_list">error list message</param>
		private void processContinuousTables(string building_type, Input input, List<Output> out_list, List<string> err_list)
		{
			Output out1 = new Output();
			List<Output> tmp_list = new List<Output>();
			
			// we need to interpolate all HVAC energy results
            List<Keyword.energy_cat_t> energy_list = createHvacEnergyList();

			// envelope table. we need to interpolate envelope's U value
			// TODO, interpolate window WWR. This function becomes a double interpolation
			tmp_list = queryTableList(building_type, input, Keyword.table_t.Envelope_TT, err_list);
			out_list[(int)Keyword.table_t.Envelope_TT] = interpolateEnergy(input, tmp_list, Keyword.keyword_t.KWT_030102, energy_list, err_list);

			// infiltration table. we need to interpolate infiltration air change rate value
			tmp_list = queryTableList(building_type, input, Keyword.table_t.Infiltration_TT, err_list);
			out_list[(int)Keyword.table_t.Infiltration_TT] = interpolateEnergy(input, tmp_list, Keyword.keyword_t.KWT_030201, energy_list, err_list);

			// window table. we need to interpolate window (WWR) rate value
			tmp_list = queryTableList(building_type, input, Keyword.table_t.Window_TT, err_list);
			out_list[(int)Keyword.table_t.Window_TT] = interpolateEnergy(input, tmp_list, Keyword.keyword_t.KWT_030401, energy_list, err_list);

			// roof table. we need to interpolate roof's U value
			tmp_list = queryTableList(building_type, input, Keyword.table_t.Roof_TT, err_list);
			out_list[(int)Keyword.table_t.Roof_TT] = interpolateEnergy(input, tmp_list, Keyword.keyword_t.KWT_030302, energy_list, err_list);

			// boilor table. we need to interpolate boilor COP value.
			tmp_list = queryTableList(building_type, input, Keyword.table_t.Heating_TT, err_list);
			// boilor COP
			out_list[(int)Keyword.table_t.Heating_TT] = interpolateEnergy(input, tmp_list, Keyword.keyword_t.KWT_050101, energy_list, err_list);

			// boilor table. we need to interpolate boilor COP value.
			tmp_list = queryTableList(building_type, input, Keyword.table_t.Cooling_TT, err_list);
			// chiller COP
			out_list[(int)Keyword.table_t.Cooling_TT] = interpolateEnergy(input, tmp_list, Keyword.keyword_t.KWT_050201, energy_list, err_list);
		}


        /// <summary>
        /// process building envelope table for continuous variable wall U and infiltration. Bilinear interpolation is called
        /// </summary>
        /// <param name="input">input object</param>
        /// <param name="out_list">energy output list</param>
        /// <param name="err_list">error mseeage list</param>
/*		private void processEnvelopeTable(Input input, List<Output> out_list, List<string> err_list)
        {
            Output out1 = new Output();
			List<List<Output>> result_list = new List<List<Output>>();

            // we need to interpolate all HVAC energy results
            List<Keyword.energy_cat_t> energy_list = createHvacEnergyList();

			// start query.
			List<Output> tmp_list = new List<Output>();

			Input new_input = new Input();

			// we have C2 fixed at lower bound, and only see C3's lower and upper bound values.
			for (int i = 0; i < (int)Keyword.keyword_t.NUM_KWT; ++i)
			{
				new_input.m_groupIndex[(int)Input.InputBound_T.Lower_IBT][i] = input.m_groupIndex[(int)Input.InputBound_T.Lower_IBT][i];
				new_input.m_groupIndex[(int)Input.InputBound_T.Upper_IBT][i] = input.m_groupIndex[(int)Input.InputBound_T.Upper_IBT][i];
			}
			new_input.m_groupIndex[(int)Input.InputBound_T.Upper_IBT][(int)Keyword.keyword_t.C2_KWT] = new_input.m_groupIndex[(int)Input.InputBound_T.Lower_IBT][(int)Keyword.keyword_t.C2_KWT];

			// query continuous variables
			tmp_list = queryTableList(new_input, Keyword.table_t.Envelope_TT, err_list);

			result_list.Add(tmp_list);

			// reset upper bound variable index
			new_input.m_groupIndex[(int)Input.InputBound_T.Upper_IBT][(int)Keyword.keyword_t.C2_KWT] = input.m_groupIndex[(int)Input.InputBound_T.Upper_IBT][(int)Keyword.keyword_t.C2_KWT];
			// we have C2 fixed at upper bound, and only see C3's lower and upper bound values.
			new_input.m_groupIndex[(int)Input.InputBound_T.Lower_IBT][(int)Keyword.keyword_t.C2_KWT] = new_input.m_groupIndex[(int)Input.InputBound_T.Upper_IBT][(int)Keyword.keyword_t.C2_KWT];
			// query continuous variables again to get another output list with dimension 2.
			tmp_list = queryTableList(new_input, Keyword.table_t.Envelope_TT, err_list);

			result_list.Add(tmp_list);

			List<Keyword.keyword_t> kw_list = new List<Keyword.keyword_t>();
			kw_list.Add(Keyword.keyword_t.C2_KWT);
			kw_list.Add(Keyword.keyword_t.C3_KWT);

			out_list[(int)Keyword.table_t.Envelope_TT] = interpolateEnvelope(input, result_list, kw_list, energy_list, err_list);
		}
		*/

        /// <summary>
        /// create energy category list for HVAC and return. The list will be used for interpolate energy based on its category.
        /// </summary>
        /// <returns></returns>
        private static List<Keyword.energy_cat_t> createHvacEnergyList()
        {
            // we need to interpolate all HVAC energy results
            List<Keyword.energy_cat_t> energy_list = new List<Keyword.energy_cat_t>();
            energy_list.Add(Keyword.energy_cat_t.Cooling_ECT);
            energy_list.Add(Keyword.energy_cat_t.CoolingTower_ECT);
            energy_list.Add(Keyword.energy_cat_t.Fan_ECT);
            energy_list.Add(Keyword.energy_cat_t.Heating_ECT);
            energy_list.Add(Keyword.energy_cat_t.Pump_ECT);

            return energy_list;
        }


        //** ---------- interpolation section ------------ **//

		/// <summary>
		/// bilinear interloplation section to interpolate envelope U value and infiltration
		/// </summary>
		/// <param name="input">input object contain input index and string</param>
		/// <param name="out_list">Energy query Output object list</param>
		/// <param name="kws">Keyword list indicating which keyword we are interpolating</param>
		/// <param name="energy_list">Energy category list indicating which energy category results we are interpolating</param>
		/// <param name="err_list">error list</param>
		/// <returns>the interpolated output object</returns>
		private static Output interpolateEnvelope(Input input, List<List<Output>> out_list, List<Keyword.keyword_t> kws, List<Keyword.energy_cat_t> energy_list, List<string> err_list)
		{
			// in default, we use lower bound result for our output.
			Output out1 = new Output();
			out1 = out_list[(int)Input.InputBound_T.Lower_IBT][(int)Input.InputBound_T.Lower_IBT];
			
			// for bilinear interpolation, we fix two continuous variable for interpolation.
            if (out_list.Count != m_bilinearNum)
            {
                err_list.Add("Dimension two for bilinear interpolation is needed. But the input variable dimension does not meet this standard.");
                return out1;
            }

            if (out_list[0].Count != (int)Input.InputBound_T.NUM_IBT || out_list[1].Count != (int)Input.InputBound_T.NUM_IBT)
            {
                err_list.Add("Energy output results should include lower and upper two object. The list number is not equal to dimension " + (int)Input.InputBound_T.NUM_IBT + ".");
                return out1;
            }

            // the nunmber of keyword should be equal to bilinear interpolation number.
            if (kws.Count != m_bilinearNum)
            {
                err_list.Add("The keyword to be interpolated number should be equal to two. Please check your input keyword number");
                return out1;
            }

            // if the keyword is direcreate variable, we don't need to do interpolation
            if (interf.BEE_enumMeasureInputType[(int)kws[0]] == interf.BEE_Input_T.Discrete_IT || interf.BEE_enumMeasureInputType[(int)kws[1]] == interf.BEE_Input_T.Discrete_IT)
            {
                err_list.Add("Can not interpolate discrete ECM variable. Please check your interpolation method.");
                return out1;
            }

            double[] X = new double[m_bilinearNum];
            double[] Y = new double[m_bilinearNum];

            double x = Convert.ToDouble(input.m_strInput[(int)kws[0]]);
			int x1_idx = input.m_groupIndex[(int)Input.InputBound_T.Lower_IBT][(int)kws[0]];
            X[0] = Convert.ToDouble(Keyword.m_inputs[(int)kws[0]][x1_idx]);
			int x2_idx = input.m_groupIndex[(int)Input.InputBound_T.Upper_IBT][(int)kws[0]];
            X[1] = Convert.ToDouble(Keyword.m_inputs[(int)kws[0]][x2_idx]);

            double y = Convert.ToDouble(input.m_strInput[(int)kws[1]]);
			int y1_idx = input.m_groupIndex[(int)Input.InputBound_T.Lower_IBT][(int)kws[1]];
			Y[0] = Convert.ToDouble(Keyword.m_inputs[(int)kws[1]][y1_idx]);
			int y2_idx = input.m_groupIndex[(int)Input.InputBound_T.Upper_IBT][(int)kws[1]];
			Y[1] = Convert.ToDouble(Keyword.m_inputs[(int)kws[1]][y2_idx]);

            // the energy result array
            double[][] Q = new double[][] 
            { 
                new double[(int)Input.InputBound_T.NUM_IBT],
                new double[(int)Input.InputBound_T.NUM_IBT]
            };

            // loop all energy type we are going to interpolate.
            for (int i = 0; i < energy_list.Count; ++i)
            {
                int energy_idx = (int)energy_list[i];

                Q[0][0] = out_list[0][(int)Input.InputBound_T.Lower_IBT].m_energy[energy_idx];
                Q[0][1] = out_list[0][(int)Input.InputBound_T.Upper_IBT].m_energy[energy_idx];
                Q[1][0] = out_list[1][(int)Input.InputBound_T.Lower_IBT].m_energy[energy_idx];
                Q[1][1] = out_list[1][(int)Input.InputBound_T.Upper_IBT].m_energy[energy_idx];

				// if we have one variable is the same, bilinear interlopation becomes linear interpolation.
				if (X[0] == X[1])
				{
					out1.m_energy[energy_idx] = linearInterpolate(y, Y[0], Y[1], Q[0][0], Q[0][1]);
				}
				else if (Y[0] == Y[1])
				{
					out1.m_energy[energy_idx] = linearInterpolate(x, X[0], X[1], Q[0][0], Q[1][0]);
				}
				else
				{
					out1.m_energy[energy_idx] = bilinearInterpolate(x, y, X, Y, Q);
				}
            }

            return out1;
		}


		/// <summary>
		/// interploate energy data based on input keyword
		/// </summary>
		/// <param name="input">the Input object</param>
		/// <param name="out_list">the Output object list, access via Input.InputBound_T</param>
		/// <param name="kw">the ECM keyword we are going to interpolate</param>
		/// <param name="energy_list">the energy category list to indicate what energy this function is going to interpolate</param>
		/// <param name="err_list">err list</param>
		/// <returns>The output object containing the interpolated energy data.</returns>
		private static Output interpolateEnergy(Input input, List<Output> out_list, Keyword.keyword_t kw, List<Keyword.energy_cat_t> energy_list, List<string> err_list)
		{
			// in default, we use lower bound result for our output.
			Output out1 = new Output();
			out1 = out_list[(int)Input.InputBound_T.Lower_IBT];

			// if the keyword is direcreate variable, we don't need to do interpolation
			if (interf.BEE_enumMeasureInputType[(int)kw] == interf.BEE_Input_T.Discrete_IT)
			{
				err_list.Add("Can not interpolate discrete ECM variable. Please check your interpolation method.");
				return out1;
			}

			if (out_list.Count != (int)Input.InputBound_T.NUM_IBT)
			{
				err_list.Add("The energy output results to be interpolated should be dimension: " + (int)Input.InputBound_T.NUM_IBT +
					". Please check the energy results list to be interpolated");
				return out1;
			}

			double x = Convert.ToDouble(input.m_strInput[(int)kw]);
			int x1_idx = input.m_groupIndex[(int)Input.InputBound_T.Lower_IBT][(int)kw];	// we store index for query table, but to get original index, we do -1
			double x1 = Convert.ToDouble(Keyword.m_inputs[(int)kw][x1_idx]);
			int x2_idx = input.m_groupIndex[(int)Input.InputBound_T.Upper_IBT][(int)kw] ;	// we store index for query table, but to get original index, we do -1
			double x2 = Convert.ToDouble(Keyword.m_inputs[(int)kw][x2_idx]);

			// loop all energy type we are going to interpolate.
			for (int i = 0; i < energy_list.Count; ++i )
			{
				int energy_idx = (int)energy_list[i];
				double y1 = out_list[(int)Input.InputBound_T.Lower_IBT].m_energy[energy_idx];
				double y2 = out_list[(int)Input.InputBound_T.Upper_IBT].m_energy[energy_idx];
				out1.m_energy[energy_idx] = linearInterpolate(x, x1, x2, y1, y2);
			}

			return out1;
		}


		/// <summary>
		/// linear interpolation, return the interpolated y value. We make it public for cost interpolation.
		/// </summary>
		/// <param name="y1">y1</param>
		/// <param name="y2">y2</param>
		/// <param name="mu">the disctance ratio defined as (x-x1)/(x1-x2)</param>
		/// <returns>the interpolated y value</returns>
		private static double linearInterpolate(double y1, double y2, double mu)
		{
			return (y1 * (1 - mu) + y2 * mu);
		}


		/// <summary>
		/// linear interpolation, return the interpolated y value
		/// </summary>
		/// <param name="x">the point's x value we are going to interpolate</param>
		/// <param name="x1">point 1, x1</param>
		/// <param name="x2">point 2, x2</param>
		/// <param name="y1">point 1, y1</param>
		/// <param name="y2">point 2, y2</param>
		/// <returns>the interpolated y value</returns>
		public static double linearInterpolate(double x, double x1, double x2, double y1, double y2)
		{
			if (x2 == x1)
				return y1;

			double mu = Math.Abs((x - x1) / (x2 - x1));
			return (y1 * (1 - mu) + y2 * mu);
		}


		/// <summary>
		/// overloading linear interpolation function
		/// </summary>
		/// <param name="x">input x point</param>
		/// <param name="X">X coordinate array X = [x1, x2]</param>
		/// <param name="Y">Y coordinate array Y = [y1, y2]</param>
		/// <returns></returns>
		private static double linearInterpolate(double x, double[] X, double[] Y)
		{
			if (X[0] == X[1])
				return Y[0];

			double mu = Math.Abs((x - X[0]) / (X[1] - X[0]));
			return (Y[0] * (1 - mu) + Y[1] * mu);
		}


		/// <summary>
		/// from wikipedia. http://en.wikipedia.org/wiki/Bilinear_interpolation
		/// </summary>
		/// <param name="x">the input x coordinate we are going to interpolate</param>
		/// <param name="y">the input y coordinate we are going to interpolate</param>
		/// <param name="X">the input x coordinate array for lower and upper boundary</param>
		/// <param name="Y">the input y coordinate array for lower and upper boundary</param>
		/// <param name="Q">the resulst array (2*2), sotred in [y][x] style</param>
		/// <returns>the interpolated result</returns>
		private static double bilinearInterpolate(double x, double y, double[] X, double[] Y, double[][] Q)
		{
			double denorm = (X[1] - X[0]) * (Y[1] - Y[0]);

			double p =	  Q[0][0] * (X[1] - x) * (Y[1] - y) / denorm
						+ Q[1][0] * (x - X[0]) * (Y[1] - y) / denorm
						+ Q[0][1] * (X[1] - x) * (y - Y[0]) / denorm
						+ Q[1][1] * (x - X[0]) * (y - Y[0]) / denorm;

			return p;
		}


		/// <summary>
		/// overloading fuction, query continuous variable based on its lower and upper bound value.
		/// use this function for specific table type and continous variable which need to query lower and upper bound.
		/// </summary>
		/// <param name="input"></param>
		/// <param name="table"></param>
		/// <param name="err_list"></param>
		/// <returns></returns>
		private List<Output> queryTableList(string building_type, Input input, Keyword.table_t table, List<string> err_list)
		{
			List<Output> outlist_temp = new List<Output>();

			// current error message number
			int err_size = err_list.Count;

			// query each table type
			for (int j = 0; j < (int)Input.InputBound_T.NUM_IBT; ++j)
			{
				// the output we would like to store to for each table
				Output out1 = new Output();

				// query original input, add if to avoid redudant input
				if (err_list.Count == err_size)
				{
					m_query.queryTable(table, building_type, input.m_groupIndex[j], out1, err_list);
				}

				outlist_temp.Add(out1);
			}

			return outlist_temp;		
		}


		/// <summary>
		/// overloading function, query all table types and add the output results into a list for return.
		/// use this function for discrete value and all table type list query
		/// </summary>
		/// <param name="input"></param>
		/// <param name="err_list"></param>
		/// <returns></returns>
		private List<Output> queryTableList(string building_type, Input input, List<string> err_list)
		{
			List<Output> outlist_temp = new List<Output>();

			// current error message number
			int err_size = err_list.Count;

			// query each table type
			for (int j = 0; j < (int)Keyword.table_t.NUM_TT; ++j)
			{
				Keyword.table_t table = (Keyword.table_t)j;

				// the output we would like to store to for each table
				Output out1 = new Output();

				// query original input, add if to avoid redudant input
				if (err_list.Count == err_size)
				{
					m_query.queryTable(table, building_type, input.m_inputIndex, out1, err_list);
				}

				outlist_temp.Add(out1);
			}

			return outlist_temp;
		}


        //** ----------- sizing building energy ---------- **//

        /// <summary>
        /// sizing energy output based on building's floor area and conditioned floor area
        /// </summary>
        /// <param name="type">building type, access by using enum interf.BEE_Building_T</param>
        /// <param name="floor_area">building floor area, access by using enum Building.BEE_BuildingFloorArea_T</param>
        /// <param name="original">original energy output list</param>
        /// <returns>the output object list</returns>
        public static Output sizingBuilding(interf.BEE_Building_T type, double[] floor_area, Output out_original)
        {
            List<Output> out_list = new List<Output>();

            // TODO: error handling when input array does not meet with standard.

            // floor area ratio compared with baseline hard coded floor areas
            double[] ratio = new double[(int)Building.BEE_BuildingFloorArea_T.NUM_BFAT];

            //Console.WriteLine("ratio: ");
            for (int i = 0; i < (int)Building.BEE_BuildingFloorArea_T.NUM_BFAT; ++i)
            {
                ratio[i] = floor_area[i] / Building.BEE_Area[(int)type][i];
                //Console.WriteLine(ratio[i] + ", ");
            }

//            for (int out_num = 0; out_num < out_original.Count; ++out_num)
//            {
                Output out1 = new Output();

                for (int i = 0; i < (int)Keyword.energy_cat_t.NUM_ECT; ++i)
                {
                    // for HVAC energy, we scale based on conditioned floor area
                    if (i <= (int)Keyword.energy_cat_t.CoolingTower_ECT)
						out1.m_energy[i] = out_original.m_energy[i] * ratio[(int)Building.BEE_BuildingFloorArea_T.Conditioned_BFAT];
                    // exterior lighitng, squre root of ratio
                    else if (i == (int)Keyword.energy_cat_t.ExtLighting_ECT)
						out1.m_energy[i] = out_original.m_energy[i] * ratio[(int)Building.BEE_BuildingFloorArea_T.Facade_BFAT];
					// TODO: refrigeration area.
//					else if (i == (int)Keyword.energy_cat_t.Refrigeration_ECT)
//						out1.m_energy[i] = original[out_num].m_energy[i] * ratio[(int)Building.BEE_BuildingFloorArea_T.Refrigeration_BFAT];
                    // for non-HVAC energy, we scale based on total floor area
                    else
						out1.m_energy[i] = out_original.m_energy[i] * ratio[(int)Building.BEE_BuildingFloorArea_T.Total_BFAT];
                }

//                out_list.Add(out1);
//            }

            return out1;
        }


		//** ---------------- economic calculation in below ----------------- **//

		/// <summary>
		/// calculate one ECM or whole building's payback period.
		/// </summary>
		/// <param name="economic">Total energy saving (+ for saving, - for not saving), and cost data (material + labor),
		/// access via Cost.BEE_EconomyResult_T
		/// <param name="err_list">error list</param>
		public static double[] payBackPeriod(double[] eco_opt, double[] cost, double[] economic, List<string> err_list)
		{
			double[] result = new double[(int)Cost.BEE_EcmCostCat_T.NUM_ECMCCT];

            // double[] results = economic;
			if (economic.Length != (int)Cost.BEE_EconomyResult_T.NUM_ERT)
			{
				err_list.Add("Error, the input economy inputs has different dimension as program required. The input should be dimension: " + (int)Cost.BEE_EconomyResult_T.NUM_ERT);
                return economic;
			}

			// energy saving over user input evaluation years.
			result[(int)Cost.BEE_EconomyResult_T.EnergySaving_ERT] = 
				calculateDiscountSum(economic[(int)Cost.BEE_EconomyResult_T.EnergySaving_ERT], eco_opt[(int)Building.BEE_BuildingEconomic_T.DiscountRate_BET], eco_opt[(int)Building.BEE_BuildingEconomic_T.EvaluationPeriod_BET]);

			// if we don't have an energy efficiency retrofit
			if (economic[(int)Cost.BEE_EconomyResult_T.EnergySaving_ERT] > 0)
				result[(int)Cost.BEE_EconomyResult_T.PayBackPeriod_ERT] = 999999;	// take forever to get your return back.
			else
			{
				// TODO, this is debug only when we don't have full cost data.
				// we give an arbitary saving ratio and back calculate cost data from energy saving.
//				if (economic[(int)Cost.BEE_EconomyResult_T.TotalCost_ERT] < 0 && economic[(int)Cost.BEE_EconomyResult_T.EnergySaving_ERT] < 0)
//					economic[(int)Cost.BEE_EconomyResult_T.TotalCost_ERT] = -economic[(int)Cost.BEE_EconomyResult_T.EnergySaving_ERT] * 0.8;

				// we must produce positive energy saving
				if (-economic[(int)Cost.BEE_EconomyResult_T.EnergySaving_ERT] - cost[(int)Cost.BEE_EcmCostCat_T.OM_ECCT] <= 0)
					result[(int)Cost.BEE_EconomyResult_T.PayBackPeriod_ERT] = 0;
				// if this is only an O&M measure which has no labor and material cost
				else if (cost[(int)Cost.BEE_EcmCostCat_T.Labor_ECCT] + cost[(int)Cost.BEE_EcmCostCat_T.Material_ECCT] == 0)
				{
					result[(int)Cost.BEE_EconomyResult_T.PayBackPeriod_ERT] = cost[(int)Cost.BEE_EcmCostCat_T.OM_ECCT] / -economic[(int)Cost.BEE_EconomyResult_T.EnergySaving_ERT];
					if (result[(int)Cost.BEE_EconomyResult_T.PayBackPeriod_ERT] > 1)	// O&M payback muse be less than one year, since O&M is performed anully.
						result[(int)Cost.BEE_EconomyResult_T.PayBackPeriod_ERT] = 999999;
				}
				else
					result[(int)Cost.BEE_EconomyResult_T.PayBackPeriod_ERT] = (cost[(int)Cost.BEE_EcmCostCat_T.Labor_ECCT] + cost[(int)Cost.BEE_EcmCostCat_T.Material_ECCT]) /
														(-economic[(int)Cost.BEE_EconomyResult_T.EnergySaving_ERT] - cost[(int)Cost.BEE_EcmCostCat_T.OM_ECCT]);
			}

			// total cost
//			result[(int)Cost.BEE_EconomyResult_T.TotalCost_ERT] = economic[(int)Cost.BEE_EconomyResult_T.TotalCost_ERT];

			// simple payback calculation
			// since energy saving is "-", we add a "-" for energy saving.
			result[(int)Cost.BEE_EconomyResult_T.TotalCost_ERT] = (cost[(int)Cost.BEE_EcmCostCat_T.Labor_ECCT] + cost[(int)Cost.BEE_EcmCostCat_T.Material_ECCT]) +
					+ calculateDiscountSum(cost[(int)Cost.BEE_EcmCostCat_T.OM_ECCT], eco_opt[(int)Building.BEE_BuildingEconomic_T.DiscountRate_BET], eco_opt[(int)Building.BEE_BuildingEconomic_T.EvaluationPeriod_BET]);

/*			// ==================== discount rate calculation ===================== //

			// we consider discount rate calculation.
			double d = eco_opt[(int)Building.BEE_BuildingEconomic_T.DiscountRate_BET];

			// this is a complicated equation derived from series sum calculation.
			double c = (-economic[(int)Cost.BEE_EconomyResult_T.EnergySaving_ERT] - cost[(int)Cost.BEE_EcmCostCat_T.OM_ECCT]) * (1 + d) / d;
			double first_invest = cost[(int)Cost.BEE_EcmCostCat_T.Material_ECCT] + cost[(int)Cost.BEE_EcmCostCat_T.Labor_ECCT];

			// for O&M cost only
			if (first_invest <= 0)
			{
				result[(int)Cost.BEE_EconomyResult_T.PayBackPeriod_ERT] = cost[(int)Cost.BEE_EcmCostCat_T.OM_ECCT] / -economic[(int)Cost.BEE_EconomyResult_T.EnergySaving_ERT];
				// for O&M cost, we must yield a number less than 1. Otherwise, O&M will have no benefit.
				if (result[(int)Cost.BEE_EconomyResult_T.PayBackPeriod_ERT] > 1)
					result[(int)Cost.BEE_EconomyResult_T.PayBackPeriod_ERT] = 999;
			}
			else if (d > 0 && result[(int)Cost.BEE_EconomyResult_T.PayBackPeriod_ERT] > 1)
			{
				// for cost with investment and labor
				result[(int)Cost.BEE_EconomyResult_T.PayBackPeriod_ERT] = Math.Log((first_invest + c) / c, 1 + d);
			}

			// total cost over evaluation period
			result[(int)Cost.BEE_EconomyResult_T.TotalCost_ERT] = eco_opt[(int)Building.BEE_BuildingEconomic_T.EvaluationPeriod_BET] * cost[(int)Cost.BEE_EcmCostCat_T.OM_ECCT]
																+ cost[(int)Cost.BEE_EcmCostCat_T.Material_ECCT] + cost[(int)Cost.BEE_EcmCostCat_T.Material_ECCT];
*/
			return result;
		}


		// ** ----------------- version 2 implementation below ------------------- ** //

		public enum UValue_T
		{
			WALL_UVT,
			ROOF_UVT,
			NUM_UVT
		};

		/// <summary>
		/// 
		/// </summary>
		/// <param name="u"></param>
		/// <param name="c_insul"></param>
		/// <param name="thickness"></param>
		/// <returns></returns>
		public static double calculateNewUValue(double u, UValue_T type, double thickness)
		{ 
			double c = 0.0;
			if (type == UValue_T.ROOF_UVT)
				c = 0.034;			// hard code conductivity of roof insulation for XPS
			else
				c = 0.04;			// hard code conductivity of wall insulation for EPS

			double temp = 1 / (1 / u + thickness / c);

			return (1 / (1 / u + thickness / c));
		}

		/// <summary>
		/// Calculate infiltration rate, based on infiluence factor
		/// </summary>
		/// <param name="infil"></param>
		/// <param name="factor"></param>
		/// <returns></returns>
		public static double calculateNewInfiltration(double infil, double factor)
		{
			return infil * factor;
		}

		/// <summary>
		/// Calculate a sum net present value based on annual saving/investment
		/// </summary>
		/// <param name="value"></param>
		/// <param name="disc"></param>
		/// <param name="year"></param>
		/// <returns></returns>
		public static double calculateDiscountSum(double value, double disc, double year)
		{
			if (year <= 1)
				return value;

			double q = 1 / (1 + disc);
			return value * (Math.Pow(q, year) - 1) / (q - 1);
		}

		/// <summary>
		/// print array function simply combining all the outputs in one line.
		/// <param name="arr">the input array</param>
		/// </summary>
		private void printArray(double[] arr)
		{
			string line = "Total: " + arr.Length + " -- ";

			for (int j = 0; j < (int)arr.Length; ++j)
			{
				Console.WriteLine(arr[j] + " ;");
			}
			Console.WriteLine(line);
		}

	}
}
